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1  Introduction 

Background 

Soil  erosion  in  the  United  States  is  usually  predicted  using  a  Revised  Universal  Soil 
Loss  Equation  (Renard  et  al.  1997).  In  the  equation,  soil  loss  is  a  function  of  six  fac¬ 
tors:  rainfall-runoff  erosivity,  soil  erodibility,  topographic  features  (steepness  and 
slope  length),  ground  and  vegetation  cover,  and  support  practice.  For  a  specific 
area,  the  cover  factor  reflects  the  effect  of  ground  and  vegetation  covers  on  the 
reduction  of  soil  loss  by  reducing  runoff  and  determines  the  dynamics  of  soil  erosion 
(Benkobi  et  al.  1994).  This  factor  is  defined  as  a  function  of  ground  cover,  canopy 
cover,  and  average  minimum  height  for  rain  dripping  down  to  the  ground  from  the 
lowest  canopy  cover  layer  (Wischmeier  and  Smith  1978).  At  Fort  Richardson, 
Alaska,  various  military  training  activities  (e.g.,  off-road  vehicular  traffic)  disturb 
ground  and  vegetation  covers,  and  cause  an  increase  of  rainfall-related  runoff  and, 
thus,  soil  erosion.  Therefore,  sampling  and  mapping  this  cover  factor  in  a  cost- 
efficient  way  becomes  critical  to  monitoring  the  dynamics  of  soil  erosion. 

Ground  cover,  vegetation  cover,  and  soil  erosion  vary  spatially,  and  their  dynamics 
are  a  complex  process.  First,  military  training  activities  often  take  place  unevenly 
in  space  and  cause  spatially  variable  disturbances  of  ground  and  vegetation  covers. 
Furthermore,  the  high  frequency  of  training  activities  in  the  same  areas  disturb 
ground  and  vegetation  covers  cumulatively  and  result  in  an  increase  in  soil  erosion 
over  time.  In  other  areas  without  such  activities,  soil  erosion  decreases  due  to 
recovery  of  the  ground  and  vegetation  covers.  For  these  reasons,  in  order  to  monitor 
and  map  dynamics  of  soil  erosion,  ground  data  should  be  collected  more  in  areas 
with  high  activities  than  in  areas  with  low  activities  (Demers  2000).  Therefore,  the 
need  is  strong  to  develop  a  cost-efficient  sample  design  method.  This  method  should 
lead  to  a  variable  sampling  distance  between  plots;  that  is,  grid  spacing  that  varies 
depending  on  the  spatial  variation  of  soil  erosion. 

On  the  other  hand,  the  supports  used  to  collect  ground  data  are  often  smaller  than 
the  desirable  units  of  mapping  or  pixels  of  images  used  because  budget  and  effort 
are  limited.  For  example,  squares  of  1  m  x  1  m  or  transect  lines  of  30  to  100  m  are 
usually  used  for  collection  of  ground  data  for  ground  and  vegetation  covers,  while 
the  pixel  size  of  Landsat  Thematic  Mapper  (TM)  imagery  used  for  mapping  is  30  m 
x  30  m.  Another  example  is  that  global  land  use  and  land  cover  is  often  mapped 


2 


ERDC/CERL  TR-06-5 


using  Moderate  Resolution  Imaging  Spectroradiometer  or  Advanced  Very  High 
Resolution  Radiometer  at  the  spatial  resolution  of  1  km  x  1  km,  while  ground  data 
are  generally  collected  at  sample  plots  that  are  smaller  than  50  m  x  50  m.  The 
inconsistency  of  support  sizes  leads  to  difficulty  and  complexity  in  developing  the 
sample  design  and  mapping  procedure.  For  these  reasons,  a  sampling  and  mapping 
method  was  developed  in  this  study. 

“Sample  design”  deals  with  determining  plot  and  sample  size,  plot  shape,  and  spa¬ 
tial  allocation  of  plots  to  collect  ground  data  given  a  cost  and  desired  precision  of 
estimates.  This  study  focused  on  determining  a  cost-efficient  sample  size  given  a 
plot  size  and  map  unit  size  for  mapping  the  cover  factor.  The  cost-efficient  sample 
size  is  defined  as  the  number  of  plots  that  minimizes  cost  given  a  desired  precision 
of  an  estimate.  Generally,  a  sample  design  based  on  random  selection  of  data 
locations  often  leads  to  a  larger  sample  and  duplication  of  information  because  some 
plots  are  too  close  together.  A  sample  design  based  on  stratified  random  selection  of 
data  locations  can  increase  cost-efficiency  by  subdividing  the  total  military  base  into 
homogeneous  strata  such  that  variances  are  minimized  within  strata  (Thompson 
1992).  In  the  traditional  sample  design  and  mapping  methods,  variables  are 
assumed  to  have  no  spatial  correlation  and  also  the  models  to  collect  ground  data 
obtained  at  smaller  supports  are  often  used  directly  to  estimate  the  values  of  the 
variables  at  coarser  spatial  resolutions.  The  methods  neglect  the  spatial  correlation 
of  variables  and  change  across  spatial  resolutions. 

Geostatistical  methods  such  as  cokriging  are  based  on  the  regionalized  variable  the¬ 
ory  that  values  of  a  random  function  are  similar  when  they  are  close  to  each  other, 
and  the  similarity  decreases  as  the  distance  between  data  locations  increases 
(Curran  and  Atkinson  1998).  This  characteristic  is  also  called  spatial  autocorrela¬ 
tion  or  spatial  variability  of  a  random  function.  Curran  (1988)  suggested  the  im¬ 
provement  in  determining  a  sample  size  using  minimized  cokriging  variance-based 
sample  design.  Atkinson  et  al.  (1992,  1994,  2000)  and  McBratney  and  Webster 
(1983)  demonstrated  this  method  and  its  applications. 

Sample  designs  for  collecting  ground  data  should  be  conducted  to  reduce  the  cost  of 
mapping  the  ground  and  vegetation  cover  factor  given  a  precision  or  to  increase  the 
accuracy  of  maps  given  a  budget.  Remote  sensing  imagery  provides  complete  cover¬ 
age  of  the  ground  surface.  Because  of  significant  correlation  of  this  cover  factor 
with  spectral  variables,  remotely  sensed  data  can  be  used  to  improve  the  estimation 
of  unknown  locations  from  the  sample  data  (Wang  et  al.  2002).  The  widely  used 
approaches  include  classification  methods  such  as  supervised  and  unsupervised 
classification  or  stratification,  and  mapping  methods  such  as  regression  estimation, 
cokriging,  and  co-simulation  in  geostatistics  (Campbell  1996;  Goovaerts  1997). 
Wang  et  al.  (2002)  compared  the  methods  for  mapping  this  cover  factor  and 
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concluded  that  cokriging  with  TM  imagery  led  to  higher  accuracy  than  unsuper¬ 
vised  stratification  and  linear  regression  by  strata. 

Cokriging  with  remotely  sensed  data  can  be  used  to  map  variables  of  interest  and 
cokriging  variance-based  sample  design  can  be  applied  to  design  sampling  strate¬ 
gies  prior  to  mapping  (Curran  and  Atkinson  1998).  Cokriging  makes  use  of  spatial 
autocorrelation  and  cross  correlation  to  explore  and  describe  spatial  variability  of 
variables,  and  further  to  optimally  estimate  local  values  from  data  sampled  else¬ 
where  and  provide  the  minimized  error  or  cokriging  variances  (Goovaerts  1997). 
The  cokriging  variances  depend  on  the  spatial  configuration  of  the  data  locations 
but  not  on  the  data  values,  and  thus  can  be  used  to  determine  the  sampling  distance 
—  grid  spacing  for  a  systematical  sample  design  to  collect  ground  data  of  a  primary 
variable  (Curran  and  Atkinson  1998). 

Because  of  the  limited  budget,  however,  the  size  of  the  plots  used  to  collect  ground 
data  is  usually  smaller  than  the  pixels  of  images  used  for  mapping.  The  inconsis¬ 
tency  of  support  sizes  means  that  calculated  cokriging  variances  do  not  really  reflect 
the  uncertainties  of  estimates.  An  alternative  is  to  use  block  cokriging  to  scale  up 
the  ground  data  at  smaller  supports  for  coarser  spatial  resolutions  and  calculate 
corresponding  block  cokriging  variances  (Deutsch  and  Journel  1998;  Goovaerts 
1997). 

The  block  cokriging  variance-based  sample  designs  are  not  locally  optimal,  however, 
because  the  error  variances  do  not  depend  on  the  data  values  and  thus  this  method 
neglects  the  heterogeneity  of  a  variable  from  area  to  area.  The  heterogeneity  re¬ 
quires  the  variable  grid  spacing;  that  is,  more  ground  measurements  need  to  be  col¬ 
lected  in  the  parts  of  a  study  area  with  higher  variation  than  in  the  parts  with  lower 
variation.  The  potential  solution  for  this  requirement  may  be  to  first  stratify  the 
study  area  into  homogeneous  strata  and,  within  each  stratum,  the  spatial  variabil¬ 
ity  of  variables  is  then  modeled  and  used  for  sample  designs.  This  method  will 
make  use  of  the  advantages  of  both  stratification  and  minimized  block  cokriging 
variance  for  sample  designs.  Integrating  stratification  and  block  cokriging  vari¬ 
ance-based  sampling  design  will  thus  increase  the  cost-efficiency  for  sampling  and 
mapping. 


Objective 


The  objective  of  this  study  was  to  develop  a  cost-efficient  method  to  sample  and  map 
the  cover  factor  when  the  size  of  supports  used  to  collect  ground  data  is  smaller 
than  the  spatial  resolution  of  map  units.  This  method  should  lead  to  a  sample 
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design  with  variable  grid  spacing  that  varies  depending  on  spatial  variation  of  the 
cover  factor. 


Approach 

The  method  was  developed  by  integrating  stratification  and  block  ordinary  co- 
kriging  with  TM  imagery.  The  TM  images  help  increase  the  cost-efficiency  of 
sampling  and  mapping  the  variable.  The  methods  were  then  applied  to  Fort 
Richardson,  where  the  cover  factor  was  sampled  and  mapped  for  monitoring  the 
dynamics  of  soil  erosion  due  to  military  training  activities. 


Scope 

This  report  analyzes  to  U.S.  Army  Integrated  Training  Area  Management  (IT AM) 
Range  and  Training  Land  Assessment  (RTLA)  program  data  from  Fort  Richardson. 
The  methodology  is  applicable  to  all  installations  involved  in  the  ITAM  RTLA  pro¬ 
gram.  While  different  installations  use  various  field  methods  to  collect  ground  data, 
the  data  extrapolation  methods  described  in  this  report  are  applicable  to  all  of  these 
installations. 


Mode  of  Technology  Transfer 

The  information  in  this  report  will  be  provided  to  U.S.  Army  ITAM  personnel 
responsible  for  RTLA  program  implementation. 

This  report  will  be  made  accessible  through  the  World  Wide  Web  (WWW)  at  URL: 
http://www.cecer.army.mil 
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2  Study  Area  and  Data  Sets 


The  study  area  consists  of  25,184  ha  (62,229  acres)  located  at  Fort  Richardson  in 
south-central  Alaska  adjacent  to  the  city  of  Anchorage  and  in  a  transitional  zone 
between  the  maritime  climatic  zone  to  the  south  and  the  interior  or  continental  cli¬ 
matic  zone  to  the  north  (Colorado  State  University  2004).  The  area  lies  in  an  allu¬ 
vial  plain  that  is  bordered  on  the  east  by  the  Chugach  Mountains  and  on  the  north, 
south,  and  west  by  waters  of  Cook  Inlet.  Figure  1  shows  a  digital  elevation  model 
(DEM)  and  slope  map  of  this  area.  Within  this  area,  elevation  decreases  from 
southeast  to  northwest  and  varies  from  16  to  5,329  m,  and  the  range  of  slope  at  the 
spatial  resolution  of  30  m  x  30  m  changes  from  0  to  82  degrees.  Average  annual 
precipitation  is  388  mm  (15.29  in.)  with  an  average  annual  snowfall  of  1,956  mm  (77 
in.).  The  area  is  dominantly  covered  by  forests  (55.3  percent),  then  scrub  lands 
(23.7  percent),  human  disturbed  lands  (13.1  percent),  barren  lands  (5.5  percent),  bog 
and  wetland  (1.6  percent),  meadow  (0.7  percent),  and  water  (0.5  percent)  (Jorgensen 
et  al.  2002).  Lowland  black  spruce  forests  are  common  and  bottomland  riparian  ar¬ 
eas  typically  have  spruce/poplar  forests  along  with  willow  and  alder  shrubs.  Mo¬ 
raines  support  white  spruce  forests.  Cottonwood/tall  bush  communities  are  common 
on  the  floodplains.  Ground  and  vegetation  covers  were  disturbed  due  to  various 
military  training  activities  from  1997  to  2001,  making  soil  erosion  a  big  concern. 

For  the  dynamics  of  soil  loss,  different  numbers  of  sample  plots  were  measured  from 
1997  to  2001,  and  the  locations  of  the  sample  plots  also  varied  from  year  to  year, 
depending  on  the  areas  in  which  military  training  activities  took  place.  In  the 
southeast  of  the  study  area,  there  were  no  or  few  plots  measured  because  of  steep¬ 
ness,  high  cost  of  sampling,  and  no  activities.  The  stratified  random  sample  design 
was  used.  That  is,  the  study  area  was  first  segmented  into  homogeneous  polygons 
based  on  vegetation  and  soil  types,  and  the  polygons  were  randomly  drawn.  Within 
each  of  the  selected  polygons,  several  transects  were  again  selected  at  random. 
Along  each  of  the  selected  transects,  a  certain  number  of  1  x  1  m2  square  plots  were 
designed  and  within  each  of  the  square  plots,  ground  cover,  canopy  cover,  and 
minimum  rain  drip  vegetation  heights  were  recorded.  The  ground  and  vegetation 
cover  percentages  of  each  plot  were  calculated.  The  soil  erosion  relevant  ground 
and  vegetation  cover  factor  was  then  derived  for  each  plot  using  the  empirical  mod¬ 
els  described  by  Wischmeier  and  Smith  (1978). 
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Digital  Elevation  Model  (DEM)  Slope 
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Figure  1.  Digital  elevation  model  (DEM)  and  slope  map  of  study  area. 


In  this  study,  a  ground  data  set  consisting  of  774  plots  measured  in  1999  was  used 
to  develop  the  method  described  above.  The  locations  and  cover  factor  values  of  the 
sample  plots  are  shown  in  Figure  2,  and  the  values  vary  from  0.006  to  0.45.  In  addi¬ 
tion,  a  scene  of  Landsat  7  ETM+  images  at  a  spatial  resolution  of  30  m  x  30  m  was 
acquired  for  the  year  1999.  These  images  were  geo-referenced  to  Universal  Trans¬ 
verse  Mercator  (UTM)  using  orthophotographs  and  affin  transformation.  Six  con¬ 
trol  points  were  selected  and  the  root  mean  square  error  (RMSE)  obtained  was 
8.98  m. 

In  addition  to  original  bands,  a  normal  difference  vegetation  index  and  other  four- 
ratio  images  were  calculated.  The  ratio  images  were  TM5/TM4,  TM7/TM4, 
(TM3+TM5)/TM4,  and  (TM3+TM7)/TM4.  The  coefficients  of  correlation  between  the 
cover  factor  and  the  images  were  then  calculated.  Figure  2  shows  the  ratio  image 
(TM3  +  TM7)/TM4.  This  image  had  the  highest  correlation  with  the  cover  factor, 
0.47166,  and  it  was  combined  with  the  set  of  ground  data  for  sample  designs  and 
mapping  of  the  cover  factor  using  block  ordinary  cokriging. 
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Sample  plots 
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Figure  2.  Locations  and  cover  factor  values  of  sample  plots,  and  ratio  image  (TM3+TM7)/TM4  that 
has  highest  correlation  with  cover  factor. 
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3  Method 


Spatial  Correlation 

Cokriging  with  TM  images  for  mapping  and  use  of  cokriging  variance  for  sample 
design  are  based  on  the  theory  of  regionalized  variables  (Goovaerts,  1997;  Curran 
and  Atkinson  1998).  In  this  study,  the  cover  factor  and  a  spectral  variable  ratio  im¬ 
age  are  denoted  with  Z,  and  Z2 ,  respectively.  Suppose  that  u  is  a  location  in  2- 
dimensional  space,  Zj(u)  and  Z2(u)  are  random  functions.  The  random  functions 

satisfy  the  intrinsic  hypothesis  (Mather  on  1971)  if,  in  addition  to 
E[Z .  (u)  -  Z(  (u  +  h)]  =  0  for  all  u  and  h  separation  vector  of  data  locations,  the  follow¬ 
ing  condition  is  met: 


y77 ,  ( h )  =  0.5  Var[Zi  (a)  -  Z(  (a  +  h )]  "  Var  "  means  var  iance  (Eq  1) 


y77  (h)  exists  and  is  a  function  of  h  only  and  not  a  function  of  u.  The  expected 
squared  difference  yzz  ( h )  that  quantifies  spatial  autocorrelation  of  the  random 

function  is  called  the  variogram.  In  practice,  the  variogram  is  finite  for  all  u  and  h, 
and  can  be  estimated  using  the  experimental  variogram: 


1  N(h) 

fZlZi  (h)  =  Y^(h)  £  (Z‘ (U" )  “  Z' (l<"  +  h))2  (Eq  2) 

where  N(h)  is  the  number  of  all  pair-wise  Euclidean  distances,  and  z;.(wa)  and 
z.  ( ua  +  h )  are  observations  of  the  variable  Z(.  at  spatial  locations  ua  and  ua  +  h ,  re¬ 
spectively.  If  the  random  function  is  not  only  intrinsically  stationary  but  also  sec¬ 
ond-order  stationary  (i.e.,  has  both  constant  mean  and  variance),  then  its  spatial 
covariance  C77  (h)  function  exists  and  has  this  relationship  with  the  variogram: 


Czz(h)  =  Cov{Zi(u),Zi(u  +h)\  "Cov"  means  covariance 


(Eq  3) 
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cZiZl(h)  =  cZiZi(0)-rZiZi(h ) 


(Eq  4) 


Moreover,  the  cover  factor  is  spatially  cross-correlated  with  the  spectral  variable, 
and  the  spatial  cross-correlation  can  be  quantified  by  a  cross-variogram  defined  as: 


7ZZ  ( h )  =  0.5Cov{Zl (h) - Zx (u  +  h), Z2 (u)  - Z2 ( u  +  h)} 


(Eq  5) 


It  is  symmetric  if  (h)  =  y^  (-/?)  and  yZ{Zi  (h)  =  yZA  (h) . 


The  cross-variogram 


can  be  estimated  using  the  experimental  cross-variogram: 


y  (h)  =  , , 

2  N(h)U 


^  N(h) 

Z  ( Z1  “  Z1  (“a  +  h^Z2  (Ua  )  -  Z2  (U a  +  k)) 


(Eq  6) 


Under  the  second-order  stationary  assumption,  there  will  also  be  a  cross-covariance 
function  Czz  (h),  and  its  relationship  with  the  cross  variogram  yzz  (h)  when 

QiZ2  (h)  ^  CZ,Z|  (h)  is. 


CZlZ,  (h)  =  Cov{Zx  (; u ),  Z2  (m  +  h)}  (Eq  7) 

yZlz2  (h)  =  0.5[CZiZ2  (0)  +  CZjZi  (0)]  -  0.5[CZ]Z2  (h)  +  (h)]  (Eq  8) 


Linear  Model  of  Coregionalization 

For  cokriging,  a  matrix  variogram  in  which  the  diagonal  entries  are  variograms  and 
the  off-diagonal  entries  are  cross-variograms  must  be  strictly  conditionally  negative 
definite.  Meeting  this  requirement  is  not  easy  for  a  specific  function.  In  practice,  a 
limited  list  of  known  valid  models  (i.e.,  nugget  effect,  spherical,  exponential,  and 
Gaussian  models),  together  with  positive  linear  combinations  of  the  models  (i.e., 
nested  models),  is  used  (Goovaerts  1997).  This  practice  does  not  work  in  the  multi¬ 
variate  case  (i.e.,  cokriging).  A  compromise  is  to  use  the  linear  coregionalization 
model  (LCM).  If  Y(u)  =  YL(u)]  is  a  vector-valued  random  function  from  L 

orthogonal  random  functions,  then  an  LCM  would  be  of  the  form  (Goovaerts  1997): 
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1= 1 


(Eq  9) 


where  each  b:  is  an  LxL  positive  definite  matrix,  and  each  function  ^(h)  is  a  valid 
(scalar)  variogram.  If  at  least  one  of  the  b , ’s  is  strictly  positive  definite,  and  all  the 
7;(h)  ’s  are  strictly  conditionally  negative  definite,  then  y (h)  is  a  strictly  condition¬ 
ally  negative  definite  matrix  function.  Thus,  the  LCM  is  a  generalization  of  the 
nested  model  in  the  scalar  case. 

Mathematically,  the  entries  b,(i,j )  of  each  matrix  b ,  correspond  to  the  sill  or  slope 
values  of  a  positive  definite  model  yt{  h)  .  For  each  /,  b,(i,j)  =  b,(j,  i) .  Thus,  it  is  suf¬ 
ficient  that  b,  (/,  i)  >  0  and  bt  (  j,  j )  >  0 ,  and  the  determinant  of  the  matrix  is  positive 
or  zero.  In  practice,  a  basic  component  [&,(/,/)]  involved  in  the  cross  variogram 

must  appear  in  the  variograms  and  the  coefficient  matrices  associated  with  the  co¬ 
regionalization  need  to  be  checked  for  positive  definiteness.  For  details  of  the  LCM, 
refer  to  Goovaerts  (1997). 


Block  Cokriging  for  Sampling  and  Mapping 

Suppose  that  a  set  of  ground  data  obtained  for  the  cover  factor  (primary  variable  Z] 
in  the  study  area),  and  the  data  of  the  ratio  image  (secondary  variable  Z2 )  is  avail¬ 
able  at  each  location.  The  ratio  image  has  a  spatial  resolution  of  30  x  30  m2,  while 
the  ground  data  set  of  cover  factors  was  collected  at  the  plots  of  1  x  1  m2.  The  cover 
factor  needs  to  be  estimated  at  all  pixels  of  30  x  30  m2.  An  alternative  for  this  pur¬ 
pose  is  to  use  a  block  cokriging  estimator  that  can  scale  up  ground  data  from  a  finer 
spatial  resolution  to  a  coarse  one  by  combining  the  ground  and  image  data  for  esti¬ 
mations  of  each  location  (Deutsch  and  Journel  1998;  Goovaerts  1997). 

Figure  3  shows  the  block  ordinary  cokriging  that  scale  up  the  ground  sample  data 
(marked  by  “x”)  from  smaller  supports  of  1  x  1  m2  to  coarser  pixels  (solid  lines)  of  30 
x  30  m2.  Within  each  pixel  the  remotely  sensed  data  are  available  and  a  pixel  can 
be  regarded  as  a  block  consisting  of  900  1  x  1  m2  supports  (dash  lines).  Each  of  the 
smaller  supports  is  first  estimated,  then  the  estimates  are  averaged  and  used  as  an 
estimate  of  the  block.  In  practice,  the  estimate  of  the  block  can  be  directly  calcu¬ 
lated  by  assigning  the  weight  of  the  block  to  each  location  of  data  at  smaller  ground 
supports  and  coarser  image  pixels  within  a  given  neighborhood  such  as  a  radius  of 
5,000  m. 


X 


Figure  3.  Block  ordinary  cokriging  with  TM  imagery  to  scale  up  ground  sample  data  (marked  by  “x”)  from 
smaller  supports  of  1x1  m2  to  a  coarser  pixel  (solid  line  square)  of  30  x  30  m2  in  which  the  remotely  sensed 
data  are  available.  The  coarser  pixel  can  be  regarded  as  a  block  consisting  of  900  smaller  1x1  m2  supports 
(dash  lines)  whose  estimates  are  obtained  and  then  averaged  as  an  estimate  of  the  block. 
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Let  v(u)  denote  a  block  centering  at  location  u\  n{(v( «)) ,  the  number  of  ground  sam¬ 
ple  data  of  the  cover  factor;  and  n2(v(u ))  ,  the  number  of  image  data  of  the  spectral 
variable,  the  ordinary  cokriging  estimator  for  the  block  v(u)  is: 

Zi(v(w))  =  X  Z  4*(V(«)»J  (Eq  10) 

i= 1  a= 1 

where  Xia(v(u))  are  the  weights  of  the  block  to  data.  This  estimator  is  unbiased  if 

«1  (!'(»))  »2  (  '’(»))  KlOOO) 

Z  A«(V(M))  =  1  and  Z  ^2ff(v(M))  =  0  ■  The  constraint  ^  \a(v(u))  =  1  implies 

a-\  a= 1  a=l 

that  not  all  of  the  weights  for  the  cover  factor  are  zero  and  at  least  one  observation 
must  be  used  for  estimation  of  each  block.  The  system  of  the  equations  to  determine 
the  weights  in  the  block  cokriging  estimator  is  obtained  by  imposing  two  con¬ 
straints:  (1)  that  the  estimator  is  unbiased,  and  (2)  that  the  variance  of  the  error  of 
estimation  is  minimal.  Using  the  solution  of  this  system,  the  block  cokriging  vari¬ 
ance  is  calculated  (Deutsch  and  Journel  1998;  Goovaerts  1997): 

<r\  0(w»  =  C- .  (v(m),  v(m))  -  Z  z  Kx  (K«))C-lZ,  >  v(«))  -  Ml  (Eq  1 1) 

i=l  a= 1 

Where  C  (v(u),  v(u))  is  the  block  covariance  of  the  cover  factor  and  is  approxi¬ 
mated  by  the  arithmetic  average  of  the  covariances  between  any  two  discretizing 
points  within  the  block  v(u).  CZ|Z  (ua,v(u))  is  the  covariance  (i  =  1)  or  cross  covari¬ 
ance  (i  =  2)  between  data  support  ua  and  the  block  v(u)  and  is  approximated  by  the 

arithmetic  average  of  the  covariances  between  the  data  support  and  the  points  dis¬ 
cretizing  the  block  v(u).  //,  is  a  Lagrange  multiplier. 

As  kriging  variance  is  for  a  single  variable,  cokriging  variance  is  used  for  sample 
designs  of  multiple  variables  (McBratney  and  Webster  1983).  According  to  Eq.  11, 
the  block  cokriging  variance  varies  depending  on  only  the  variogram  and  cross 
variogram,  spatial  configuration  of  data  locations  in  relation  to  the  block  to  be  esti¬ 
mated,  and  on  the  distances  among  the  data  locations.  It  does  not  depend  on  the 
data  values.  If  (1)  the  variogram  and  cross  variogram  are  known;  (2)  a  systematical 
and  square  grid  pattern  is  used;  and  (3)  the  numbers  of  ground  and  image  data  used 
for  estimation  of  each  block  are  fixed,  changing  grid  spacings,  then  sampling  dis¬ 
tances  of  ground  and  image  data,  respectively,  will  lead  to  different  block  cokriging 
variances  (Atkinson  et  al.  1992;  McBratney  and  Webster  1983).  Then,  the  maxi¬ 
mum  block  cokriging  variance  can  be  used  as  the  measure  of  quality  of  the  sample 
design.  From  the  diagram  of  maximum  block  cokriging  variance  against  grid  spac¬ 
ing,  a  sampling  distance  for  the  desirable  precision  requirement  can  be  determined. 
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The  block  cokriging  variance-based  sample  design  will  lead  to  a  systematic  sample 
that  ensures  that  neighboring  sample  points  are  as  far  apart  as  possible  for  a  fixed 
area.  Thus,  this  method  minimizes  the  duplication  of  information  by  maximizing 
the  distance  between  sample  points  given  a  desired  precision,  and  it  is  more  cost- 
effective  than  the  sample  design  based  on  random  selection  of  locations  (Atkinson  et 
al.  1992,  1994,  and  2000;  Curran  and  Williamson  1986). 

Cokriging  variance-based  sample  design  for  multiple  variables  is  more  complicated 
than  the  kriging  variance-based  method  for  a  single  variable  (McBratney  and  Web¬ 
ster  1983).  The  former  requires  models  of  spatial  cross-correlation  between  vari¬ 
ables,  in  addition  to  spatial  autocorrelation  functions,  and  deals  with  different 
sampling  intensities  of  the  primary  and  secondary  (spectral)  variable  and  the  ratio 
between  them.  In  the  univariate  case,  the  maximum  kriging  variance  occurs  at  the 
locations  that  are  most  away  from  the  data  locations,  and  is  located  at  the  center  of 
spatial  configuration  of  observations  given  a  spatial  pattern  such  as  square  grid.  In 
the  multivariate  case,  however,  the  position  of  the  maximum  cokriging  variance 
changes  depending  on  the  form  of  variograms  and  on  the  strength  of  cross  correla¬ 
tion,  the  sampling  interval,  and  the  ratio  of  the  numbers  of  data  locations  for  the 
two  variables  inside  the  search  neighborhood  that  matters.  Wang  et  al.  (2005)  sug¬ 
gested  that  the  average  of  cokriging  variances  could  be  applied  as  the  measure  of 
uncertainty  for  sample  design  in  order  to  avoid  searching  for  the  maximum 
cokriging  variance.  In  addition,  the  authors  also  found  that  using  neighboring 
pixels  for  image  data  to  map  the  cover  factor  led  to  more  reliable  estimates  than 
using  the  pixels  that  were  drawn  by  a  systematic  selection. 
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4  Stratified  Sample  Designs 

In  this  study,  the  average  of  cokriging  variances  within  the  study  area  was  calcu¬ 
lated  and  plotted  against  grid  spacing  given  a  square  grid.  With  a  desirable  preci¬ 
sion,  a  sampling  distance  was  then  determined  from  the  developed  diagrams. 
Moreover,  the  neighboring  pixels  around  a  location  to  be  estimated  were  directly 
used  instead  of  the  pixels  systematically  selected  for  the  image  data  used  for  calcu¬ 
lating  cokriging  variance.  In  addition,  three  sampling  intensities  were  compared 
with  ratios  of  ground-to-image  data  from  1:1  to  1:4  and  1:9  to  explore  the  effect  of 
the  image  data  to  reduce  cost  or  increase  map  accuracy  for  sampling  and  mapping 
the  cover  factor.  For  estimation  of  each  location,  12  ground  data  were  used.  The 
number  of  image  data  for  three  ratios  varied  from  12  to  48  and  108. 

Using  remotely  sensed  data,  a  study  area  usually  can  be  stratified  into  homogene¬ 
ous  subareas  or  strata.  Within  each  of  the  subareas  or  strata,  the  spatial  variability 
of  cover  factor  is  similar,  and  a  sample  design  can  be  made.  This  will  lead  to  vari¬ 
able  grid  spacing  that  is  optimal  within  the  corresponding  subarea  or  stratum.  In 
this  study,  the  southeast  part  of  the  area  has  a  slope  greater  than  45  degrees,  and 
no  or  few  military  training  activities  took  place.  The  off-road  vehicular  traffic 
mainly  occurred  in  the  southwest,  northwest,  and  northeast  flat  areas.  Instead  of 
using  remotely  sensed  data,  therefore,  this  study  area  was  directly  divided  into 
training  and  steep  subareas,  and  a  sample  design  was  then  carried  out  within  each 
of  them. 

By  combining  the  ground  data  and  ratio  image  mentioned  above,  an  estimation  map 
of  the  cover  factor  was  generated  with  a  sequential  Gaussian  collocated  co-simula¬ 
tion  (Anderson  et  al.  2005;  Wang  et  al.  2002).  In  the  co-simulation,  500  realizations 
of  the  cover  factor  map  were  generated,  and  a  sample  average  map  was  calculated 
as  the  estimation  of  the  cover  factor  surface.  The  map  was  then  used  as  a  reference 
surface,  based  on  which  a  stratified  sample  design  and  a  sample  design  without 
stratification  were  constructed  using  the  block  cokriging  variance-based  method  de¬ 
scribed  above.  Using  these  two  samples  together  with  the  original  sample,  the  cover 
factor  was  mapped  and  the  resulting  maps  were  compared  with  this  reference  map. 
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5  Results 


The  experimental  variograms  y  cp(  I  h  |)  and  y  m  (|  h  |)  of  the  cover  factor  and  ratio 


image  (TM3  +  TM7)  /  TM4,  and  experimental  cross-variogram  y  ^  (\h\)  between 

them  were  calculated  for  the  entire  area  using  VARIOWIN  software  (Pannatier 
1996).  These  variograms  were  modeled  using  the  spherical  models: 


y  (|  h  |)  =  0.0055  +  0.0069[1. 5—^ — 0.5(^^)3]  0  <|  /?  |<  12000  (Eq  12) 

/  cfvi  !2000  12000 


A  iii  |  i  | 

7  (|  h  |)  =  140.0  + 180.0[1.5  — — - 0.5(^^)3]  0  <|  /z  |<  12000  (Eq  13) 

/  tmvi  12000  12000 


/\  |  i  |  |  i  | 

y  (|  h  |)  =  0.272  +  0.792[1.5  — — - 0.5(^^)3]  0  <|  /?  |<  12000  (Eq  14) 

/  iCF,TMy  !2000  12000 


where  |  h  |  is  the  magnitude  of  separation  vector  of  data  locations.  Because  the  co¬ 
efficient  matrices  in  the  linear  coregionalization  model  must  be  positive  definite,  it 
may  be  necessary  to  adjust  the  entries  in  one  or  more  of  the  matrices.  In  that  case, 
the  resulting  variograms  may  not  fit  the  experimental  variograms  as  well  as  the 
original  fits.  Figure  4  presents  these  experimental  and  modeled  variograms  and 
cross-variograms.  The  semivariances  increased  over  the  separate  distance  given  a 
direction,  and  reached  their  maximum  values  up  to  12,000  m,  a  common  range  of 
spatial  correlation. 

With  the  models,  block  cokriging  variance  of  each  pixel  was  calculated  when  differ¬ 
ent  grid  spacings  from  120  m  to  2,580  m  (corresponding  sample  size  varied  from 
17,489  to  38  plots)  and  three  ratios  (mentioned  above)  of  ground-to-image  data  were 
used.  Figure  5  presents  the  change  in  average  value  of  the  block  cokriging  vari¬ 
ances  versus  grid  spacing  for  sampling  the  cover  factor.  For  the  same  ratio  of 
ground-to-image  data,  the  average  of  block  cokriging  variances  increased  as  the 
sampling  distance  increased. 
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Figure  4.  Experimental  and  modeled  variograms  of  (a)  cover  factor,  (b)  TM  ratio  image 
(TM3+TM7)/TM4,  and  (c)  cross  variogram  between  the  cover  factor  and  ratio  image.  y{\  h  |)  is  the 

semivariance  and  |h|  is  the  separation  distance  of  data  locations  with  the  unit  of  100  m. 


Figure  5.  Average  of  cokriging  variances  versus  grid  spacing  for  sampling 
ground  data  with  three  different  ratios  of  ground  data  to  image  data. 
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For  the  same  sampling  distance,  the  average  of  block  cokriging  variances  in  Figure 
5  was  smaller  for  the  ground-to-image  data  ratio  of  1:9,  that  is,  where  12  ground 
data  and  108  image  data  were  used,  than  for  the  ground-to-image  data  ratio  1:4, 
where  12  ground  data  and  48  image  data  were  used.  However,  the  difference  be¬ 
came  smaller  as  the  grid  spacing  increased.  For  the  ground-to-image  data  ratio  1:1 
where  12  ground  data  and  12  image  data  were  used,  the  average  of  block  cokriging 
variances  was  larger  than  those  from  both  ratios  1:4  and  1:9  when  the  grid  spacing 
was  less  than  720  m,  and  from  the  grid  spacings  of  720  m  to  1,200  m,  smaller  than 
that  from  the  ratio  1:4  but  still  larger  than  that  from  the  ratio  1:9,  and,  after  that, 
smaller  than  those  from  both  ratios  1:4  and  1:9.  The  reasons  for  the  phenomena 
above  may  include  (1)  as  a  greater  grid  spacing  was  used  to  sample  the  cover  factor, 
a  larger  block  cokriging  variance  was  obtained  from  the  spatial  configuration  of 
data  locations,  and  the  effect  of  image  data  of  neighboring  locations  in  the  reduction 
of  block  cokriging  variance  relatively  decreased;  and  (2)  as  the  separation  distance 
of  data  locations  increased,  the  spatial  cross  correlation  between  the  cover  factor 
and  spectral  variable  decreased,  that  is,  spatial  variability  increased,  which  implied 
that  using  more  image  data  than  needed  might  lead  to  more  noise. 

When  the  study  area  was  divided  into  two  strata  (steep  and  training  subareas  in 
Figure  6),  the  experimental  variograms  and  cross-variograms  were  calculated  and 
modeled  for  each  stratum.  These  variograms  include: 

Within  the  steep  subarea: 


/CF(\  h\)  =  0.000115  +  0.0004  1[1.5J^L_o.5(JAL)3]  0  <|  h  |<  6020  (Eq  16) 


y  (\h\)  =  1 15.176  +  136.792[1.5-^ — 0.5(^^)3]  0 <\h\<  6020  (Eq  17) 

'  ™  6020  6020 


v  (|  h  |)  =  0. 105  +  0.2 13[1 .5  — — 0.5(^^)3]  0<|A|<6020  (Eq  18) 

/  iCFjuy  6020  6020 


Within  the  training  subarea: 


y  (|  h  |)  =  0.0057  +  0.008 1[1 .5 — 0.5(^^)3]  0<|/*|<  12040  (Eq  19) 

/  cfvi  !2040  12040 
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A  111  111 

y  (|  A  |)  =  128.1  +  146.993[1.5— - 0.5(^^)3]  0  <|  /?  |<  12040  (Eq20) 

/  tmk'  12040  12040 


A  111  111 

y  (\h\)  =  0.294  +  0.700[1 .5  — — - — 0.5(^^)3]  0  <|  /?  |<  12040  (Eq21) 

/  ( CFJM )VI  12040  12040 
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Figure  6.  Stratification  of  study  area  into  steep  and  training  subareas. 
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The  parameters  of  the  variogram  models  differed  between  the  steep  area  to  the 
training  area.  For  example,  the  range  of  spatial  correlation  is  6,020  m  and 
12,040  m  for  the  steep  and  training  subareas,  respectively.  Figure  7  shows  the  av¬ 
erages  of  block  cokriging  variances  versus  grid  spacing  for  sampling  the  cover  factor 
with  these  two  strata.  The  ratio  1:4  of  ground-to-image  data  was  used  to  save  com¬ 
putational  time.  For  comparison,  the  global  curve  without  stratification  was  also 
presented  in  Figure  7.  As  the  grid  spacing  increased,  the  average  of  block  cokriging 
variances  increased  slowly  for  the  training  subarea  and  the  entire  area,  but  rapidly 
for  the  steep  subarea. 

Given  a  precision-average  of  block  cokriging  variance  equal  to  0.04,  Figure  7  was 
used  to  calculate  the  grid  spacing  for  the  entire  area,  training  subarea,  and  steep 
subarea,  respectively.  The  grid  spacing  obtained  was  295  m  when  the  sample  de¬ 
sign  was  carried  out  using  block  cokriging  variance-based  sampling  without  stratifi¬ 
cation,  322  m  for  the  steep  subarea,  and  765  m  for  the  training  subarea  using  block 
cokriging  variance-based  sampling  with  two  strata.  The  corresponding  sample  size 
was  2,890  without  stratification  and  938  with  stratification.  The  number  of  the 
original  sample  plots  was  747.  Figure  8  presents  these  three  samples.  Using  the 
reference  map  of  cover  factor  by  co- simulation,  the  effects  of  these  three  samples  to 
map  the  cover  factor  were  compared. 


Figure  7.  Average  of  cokriging  variances  versus  grid  spacing  for  sampling  ground  data  with  a 
ratio  1:4  of  ground-to-image  data  when  the  study  area  was  divided  into  steep  and  training 
subareas. 
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Block  cokriging  variance  based 
sample  design  without  stratification 


Block  cokriging  variance  baser 
sample  design  with  two  strata 


Original  sample  design 


Legend 
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Figure  8.  Block  cokriging-based  sample  designs  without  stratification  and  with  two  strata, 
respectively,  and  the  original  sample  design. 
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Figure  9  shows  the  estimation  maps  of  cover  factors  from  the  three  sample  designs 
just  described  using  the  block  ordinary  cokriging  estimator  with  the  ratio  image 
(TM3  +  TM7)  /  TM4  and  their  comparison  with  the  reference  map.  The  estimates  of 
the  map  by  the  original  samples  were  obviously  smoothed.  Two  maps  by  the  block 
cokriging  variance-based  sample  designs  with  and  without  stratification  were  simi¬ 
lar  to  each  other  and  also  appeared  close  to  the  reference  map  in  terms  of  spatial 
distribution  and  pattern  of  the  values. 


Cover  factor  (CF)  reference  map  CF  map  by  original  sample  design 


lillliL 

‘V-  >1 

‘Mf' 

HR 

r*V\ 

H 

CF  map  by  cokriging  variance 
based  sampling  without  stratification 


CF  map  by  cokriging  variance 
based  sampling  with  strata 


Legend 


Cover  factor 


0  7,500  15,000  Meters 


|  0.0  -  0.01 
I  0.01  -  0.04 
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|  0.08  -  0.12 
I  0.12  -  0.45 


Figure  9.  Estimation  maps  of  cover  factor  (CF)  using  three  sample  designs:  original  sample 
design,  block  cokriging  variance-based  sampling  without  stratification  and  with  two  strata 
(steep  and  training  subareas),  and  their  comparison  with  the  reference  map. 
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In  Figure  10,  the  three  sample  designs  were  further  compared  by  calculating  rela¬ 
tive  root  mean  square  errors  (RMSE  *  100  /  sample  mean)  of  these  estimation  maps 
from  the  reference  map  using  an  independent  sample  of  384  randomly  drawn  plots. 
The  original  sample  had  the  largest  relative  RMSE,  144.3  percent,  but  had  the 
smallest  sample  size,  747  plots.  The  error  of  the  sample  design  by  the  block  cok- 
riging  variance-based  sampling  with  stratification  was  slightly  larger,  73.2  percent, 
than  that  without  stratification,  62.3  percent,  while  the  sample  size  of  the  latter 
was  three  times  larger  than  that  of  the  former.  Thus,  the  cost-efficiency  per  sample 
plot  significantly  differed  between  these  sample  designs. 

Because  the  largest  relative  RMSE  was  144.3,  200  -  relative  RMSE  was  defined  as 
the  total  efficiency  of  sampling  and  mapping  and  (200  -  relative  RMSE)/sample  size 
as  the  unit  cost-efficiency.  The  block  cokriging  variance-based  sample  design  with 
stratification  had  the  highest  unit  cost-efficiency,  followed  by  the  original  sample 
design  and  the  block  cokriging  variance-based  sample  design  without  stratification. 
The  reason  that  the  original  sample  design  had  a  higher  cost-efficiency  than  the 
block  cokriging  variance-based  sample  design  without  stratification  was  that  the 
former  had  a  different  precision  requirement  and  a  significantly  smaller  sample  size 
than  the  latter.  If  the  unit  cost-efficiency  of  the  original  sample  was  regarded  as  1, 
the  unit  cost-efficiencies  by  the  block  cokriging  variance-based  sample  designs  with 
and  without  stratification  were  1.81  and  0.64,  respectively.  That  is,  the  stratifica¬ 
tion  increased  the  unit  cost-efficiency  by  three  times  when  compared  with  no  strati¬ 
fication. 
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Figure  10.  Comparison  of  three  sample  designs  for  mapping  the  cover  factor  based  on  (a)  relative 
root  mean  square  error  (RMSE)  -  (RMSE  *  100  /  mean)  and  (b)  cost-efficiency  per  plot  defined  as  (200 
-  relative  RMSE0  /  sample  size). 
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6  Conclusions  and  Discussion 

When  natural  resources,  and  environmental  and  ecological  systems  are  modeled 
and  mapped  with  the  aid  of  remotely  sensed  data,  cost-efficient  sample  designs  for 
collecting  ground  data  and  obtaining  accurate  maps  of  variables  are  widely  re¬ 
quired.  Because  of  limited  budgets,  however,  very  often  the  supports  used  to  collect 
ground  data  are  smaller  than  the  map  units  and  pixels  of  images  used.  This  incon¬ 
sistency  leads  to  difficulty  and  complexity  in  developing  the  procedure  of  sample  de¬ 
sign  and  mapping.  For  these  purposes,  a  sampling  and  mapping  method  was  devel¬ 
oped  by  integrating  stratification  and  an  up-scaling  method  —  block  cokriging 
variance-based  sample  design  with  TM  imagery.  This  method  can  lead  to  sample 
designs  with  variable  grid  spacing  and  can  be  used  to  significantly  increase  the  unit 
cost-efficiency  of  sample  data  in  mapping.  This  method  is  especially  effective  in 
cases  where  a  study  area  obviously  contains  heterogeneous  strata.  The  method  was 
validated  in  this  study  to  sample  and  map  a  ground  and  vegetation  cover  factor  for 
a  monitoring  system  of  soil  erosion  by  comparing  different  sample  designs  with  and 
without  stratification. 

The  block  cokriging  variance-based  sample  design  with  TM  imagery  is  based  on  the 
variograms  and  cross-variograms  that  provide  the  link  between  aspatial  statistics  of 
a  variable  and  geographic  space,  and  the  link  between  a  mapped  variable  and  a 
spectral  variable  geographically.  The  block  cokriging  with  TM  imagery  scales  up 
not  only  the  ground  sample  data  but  also  the  uncertainties  of  spatial  configuration 
of  locations  of  ground  and  image  data  used  and  variogram  models  from  smaller  sup¬ 
ports  to  larger  pixels  or  blocks.  For  the  same  grid  spacing  of  ground  data,  the  aver¬ 
age  of  block  cokriging  variances  generally  decreases  as  the  number  of  image  data 
used  increases.  The  decrease  becomes  slighter  and  even  disappears  with  a  further 
increased  number  of  image  data  because  the  uncertainties  and  noise  from  the  image 
data  may  cancel  out  the  effect  of  spatial  cross-correlation  between  the  primary  and 
secondary  variables.  In  addition,  using  too  many  image  data  will  lead  to  a  large  in¬ 
crease  in  computational  time  and  may  also  cause  the  cokriging  equation  system  to 
become  unstable  (Goovaerts  1997). 

The  disadvantage  of  using  block  cokriging  variance  to  determine  the  grid  spacing  of 
a  sample  is  that  the  block  cokriging  variances  do  not  reflect  the  uncertainties  from 
variation  of  the  ground  and  image  data  used.  Thus,  the  grid  spacing  obtained  is  not 
locally  optimal.  However,  the  estimation  of  a  population  mean  will  be  most  precise 
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if  the  population  can  be  partitioned  into  strata  in  such  a  way  that  within  each  stra¬ 
tum,  the  units  are  as  similar  as  possible.  That  is,  stratified  sampling  addresses 
homogeneity  of  units  within  stratum.  Thus,  integrating  stratification  and  cokriging 
variance-based  sampling  theoretically  can  lead  to  the  most  cost-efficient  sample  de¬ 
sign  with  variable  grid  spacing.  This  integration  improves  the  accuracy  and  cost 
efficiency  of  sampling  and  mapping  and  is  applicable  for  sample  designs  of  the  vari¬ 
ables  for  which  spatial  variation  varies  from  area  to  area.  This  outcome  was  veri¬ 
fied  by  the  results  of  this  case  study  in  which  the  unit  cost-efficiency  of  the  sample 
using  the  block  cokriging  variance-based  sample  design  with  stratification  was,  re¬ 
spectively,  1.8  times  higher  than  that  of  the  original  sample  using  traditionally 
stratified  sampling  and  3  times  higher  than  that  of  the  sample  using  the  block  cok¬ 
riging  variance-based  sampling  without  stratification. 
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